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Representing massless Dirac fermions on a spatial lattice poses a potential challenge known as 
the Fermion Doubling problem. Addition of a quadratic term to the Dirac Hamiltonian circum¬ 
vents this problem. We show that the modified Hamiltonian with the additional term results in 
a very small Hamiltonian matrix when discretized on a real space square lattice. The resulting 
Hamiltonian matrix is considerably more efficient for numerical simulations without sacrificing on 
accuracy and is several orders of magnitude faster than the atomistic tight binding model. Using 
this Hamiltonian and the Non-Equilibrium Green’s Function (NEGF) formalism, we show several 
transport phenomena in graphene, such as magnetic focusing, chiral tunneling in the ballistic limit 
and conductivity in the diffusive limit in micron sized graphene devices. The modified Hamiltonian 
can be used for any system with massless Dirac fermions such as Topological Insulators, opening up 
a simulation domain that is not readily accessible otherwise. 


Two-dimensional materials have attracted consider¬ 
able recent attention for their superior electronic charac¬ 
teristics and potential for electronic, opto-electronic and 
spintronic applicationJII^. Among them, graphene, topo¬ 
logical insulators and transition metal dichalcogenides 
stand out in particular. Experimental progress on such 
materials are advancing rapidl}!^ and modeling carrier 
transport to capture the new physics of these materials 
has become quite crucial. NEGF based numerical calcu¬ 
lation is widely used nowadays to accurately model nano¬ 
scale materialJ^. The formalism is extremely powerful in 
solving any quantum transport problem accurately in¬ 
cluding sophisticated contact-channel effects and various 
forms of scatterings within a self-energy correction to 
the Hamiltonian. Despite these considerable strengths, 
NEGF has so far been considered mainly as a ballistic 
quantum transport simulation platform, since it becomes 
computationally prohibitive to model diffusive systems. 
This is especially true for experimentally relevant de¬ 
vice dimensions which are often in the hundreds of nano¬ 
meters to jam regime. 

In materials such as graphene and topological insu¬ 
lators etc., electrons behave as massless Dirac fermions 
described by the Dirac Hamiltonian. Quantum trans¬ 
port simulations using the NEGF formalism require a real 
space representation of the Hamiltonian matrix to fully 
describe the channel material. While the tight-binding 
representation is valid for most practical purposes, it is 
computationally expensive. To expedite the calculation, 
a discretized version of the Dirac Hamiltonian can be 
used. It has however been shown that representing the 
Dirac Fermions on a spatial lattice poses a problem com¬ 


monly known as the Fermion Doublin^^ problem, where 
the numerical discretization of the Hamiltonian creates 
additional branches within the Brillouin zone (Fig. [^. 
One way to solve this problem is to add a quadratic term 
with the HamiltoniarP^. 

In this letter, we show that the additional term not 
only solves the Fermion Doubling problem but also has a 
profound implication for numerical simulations. We show 
that the modified Dirac Hamiltonian results in a spatial 
Hamiltonian matrix which is several orders of magnitude 
smaller than the atomistic tight binding HamiltoniarP 
while preserving the same level of accuracy over the rele¬ 
vant energy range. As a result, bandstructure and trans¬ 
port simulations using the modified Hamiltonian are sev¬ 
eral orders of magnitude faster than the tight binding, al¬ 
lowing us to carry out several large scale simulations such 
as angle dependent transmission in 1/im wide graphene 
(Fig. [^, magnetic electron focusing in a 0.4/imx0.4/im 
multi-electrode graphene (Fig. and diffusive transport 
in a 1/im x 0.5/im graphene device (Fig. . Our results 
show very good agreement with recently reported exper¬ 
imental resultJ^. Other methods such as tight-binding 
make it almost impossible to simulate devices at such 
dimension. 

The modified effective k.p Hamiltonian for graphene at 
low energy is, 

H{k) = hvF [k^a:c + kyay + p{kl + (1) 

where, up is the Fermi velocity, k = k^x + kyp is the 
wave vector, a’s are Pauli matrices representing the pseu¬ 
dospins, and fi is the reduced Planck’s constant. The last 
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FIG. 1. The first Brillouin zone of graphene band structure 
calculated using the discretized k.p Hamiltonian in Eq. (2) 
for (a) a = 0 and (b) a = 1.165. (c) Comparison of band 
structures along ky = 0 for different a and a. The lattice 
parameter a is in A. With a — 1.165, even a 2nm grid spacing 
results in a linear band structure for |Emax| ~ 0.6eV. 


term, I3{k^ + ky)az serves two purposes: (i) it circum¬ 
vents the well known Fermion doubling problenP and 
(ii) it allows us to generate a computationally efficient 
Hamiltonian using a course grid with reasonable accu¬ 
racy as shown below. The k-space Hamiltonian in Eq. 
is transformed to a real-space Hamiltonian by replacing 
kx with differential operator —idjdx^ k^ with —d'^jdx^ 
and so on. The differential operators are then discretized 
in a square lattice using the finite difference method to 
obtain, 

U = Y^ cjeci + + H.C.) 

+ X/ + H.C.j 
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where e = —AfiVYCccTz/a^ G = [idyl2a + acTz/a]^ ty = 
TiVY[—idx/2a-\-acFz/a]^ a is the grid spacing and a = 
/?/a. 

When (a = 0, Eqs. 0 and 0 reduce to the unmodified 


k.p Hamiltonian. The eigen-energy calculated from Eq. 
1^ becomes a sine function of wavevector k^ and three 
extra Dirac cones appear inside the first Brillouin zone 
as shown in Fig. f^a). This is known as the Fermion 
doubling problenJ^ The last term in Eq. [I] gets rid of 
this problem by opening bandgaps for each of these Dirac 
cones. The resulting band structure contains only one 
Dirac cone as shown in Fig. [^b). 

The effect of the extra term in Eq. 0 is more clearly 
illustrated in Fig. I^c). When a = 0 and a = 20A^, the 
band structure along = 0 is a sine function and an ex¬ 
tra Dirac cone appears at the zone boundary. This Dirac 
cone is removed when a = 1.165 and the band struc¬ 
ture closely follows the ideal band structure of graphene. 
The band structure calculated using a = 1.165 with 
a = 20A is accurate over a larger window than that is cal¬ 
culated using Of = 0 and a = lOA. Thus, the a parameter 
not only circumvents the Fermion Doubling problem but 
also enables us to use coarser lattice grid without losing 
the accuracy. With a = 20 A and a = 1.165, the Hamil¬ 
tonian size given by Eq. 0 for a 1/rm x 1/am graphene 
sheet is 5E5 x 5E5 compared to ^ 38E6 x 38E6 in atom¬ 
istic tight binding model. The band structure calculated 
using this discretized Hamiltonian closely resembles the 
ideal linear band structure within \E\ ^ 0.6eV. To ob¬ 
tain the same level of accuracy, the recently proposed 
scaled graphene modeP^ requires a Hamiltonian of size 
^ 9E6 X 9E6. 

The accuracy of the proposed Hamiltonian is demon¬ 
strated by the band structure of a 200nm graphene rib¬ 
bon shown in Fig The bands calculated using the one 
Pz orbital tight-binding model (in black) and the modi¬ 
fied Dirac Hamiltonian (in blue) are in good agreement 
in the low energy limit. For the tight binding model, the 
calculation takes about 1 hour and 30 minutes. In com¬ 
parison, the modified Hamiltonian model takes about 2 
seconds using the same number of cores. This shows that 
our proposed model is three orders of magnitude faster 
with excellent accuracy for band structure calculation. 
With the modified Dirac Hamiltonian, we employ the 
standard Recursive Green’s Function (RGF) algorithuJ^ 
and employ it for large scale simulations, in both ballistic 
and diffusive regimes. 

Fig. ^ shows the conductance calculated using this 
model (with a = 1.16 and a = 16A^) for a l/rm wide 
graphene sheet and compared with the conductance from 
the linear graphene E — K relationship. We see a good 
agreement until E^ = 0.6eV. In Fig. [^, we show 
the conductance of a electrostatically doped split-gated 
graphene pn junction, showi ng ex cellent agreement with 
the exact analytical solutiorP^^. That means the an¬ 
gle dependent transmission, an important property of 
graphene originating from its chiral nature^, is undis¬ 
torted in the modified Hamiltonian. It can be shown 
that the pseudospins of the modified Hamiltonian are 

of the form, ^ = ^1 exp*^//(/5,/cp)^ where 

f{l3, k-p) = Pk-p + y^l + P‘^k^, independent of the angle 
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FIG. 2. (a) Band structure of 200nm armchair graphene 

nanoribbon calculated using the tight binding model (in 
black) and the k.p model (in blue) showing good agreement. 
The yellow line is the linear approximation of the graphene 
band structure, (b) Conductance of a 1/im wide graphene 
sheet from NEGF simulation with the modified Hamiltonian 
along with that from linear E — K, (c) Conductance of a l/xm 
wide graphene pn junction in excellent agreement with exact 
analytical solutions, showing the model’s ability to capture 
angle dependent chiral transmission. 


{0 = tdin~^ky/kx) and therefore does not distort the chi¬ 
ral properties. 

Transverse magnetic field (TMF)s have been used in 
the past to study various transport phenomena and 
characterize surfaces and interfaceJ^. In a recent 
experimenfPl, a TMF is used to focus electrons in a mono- 
layer graphene device. The device geometry and biasing 
scheme used in our simulation are shown in Fig. ia). 
Electrons are injected from contact with a current source 
a and collected at contact b. In presence of a magnetic 
field 5, an electron follows a circular path inside the 
graphene channel with cyclotron radius Tc and is directly 
focused from contact a to contact c when 2rc = I/, where 
L is the distance between contacts a and c. This is the 
resonance condition where the voltmeter registers a large 
voltage. Thus, the magnetic field required for the focus¬ 
ing is 


B 



(3) 


where n is the density of electrons in the graphene chan¬ 
nel and m is an integer. For m > 1 the electron reaches 
contact c after skipping along the edges through multiple 
specular reflections. 

NEGF simulations for the multi-terminal device are 
shown in Fig. The contacts are modeled using the 
self-energies of semi-infinite graphene ribbons. The cur¬ 
rent at contact a is calculated using the multi-terminal 
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FIG. 3. Large scale ballistic simulation to capture mag¬ 
netic focusing in graphene, (a) Device geometry and biasing 
scheme. Electrons are injected from contact a and collected at 
contact b. The vertical magnetic field B forces the electrons 
to follow a circular path. Eor some specific magnetic fields, 
electrons are focused to contact c, resulting in a large voltage 
registered in the voltmeter, (b) The various matrix and trans¬ 
mission components that enter the NEGE simulation of the 
device, (c) The resistance, = E// as a function of carrier 
density and magnetic held showing resonances. The dashed 
lines were calculated using Eq. The channel dimensions 
are 400nm x 400nm. 


Landauer-Biittiker formalisnP^, 

= I E / dET^E) [/(/x„) - (4) 

where / is the Fermi function, /i is the electro-chemical 
potential of the contact and is the total transmission 
between contact a and contact The transmission is 
calculated using the Fisher-Lee formulsP, 

T^piE) = Tr{r„G^0T^Gfj (5) 

where is the retarded Green’s function and = 
i{'^a,i3 ~ ^ 1 / 3 ) is the broadening from contacts P) 
with being the corresponding energy dependent self¬ 
energy matrices. For computational efficiency, the re¬ 
tarded Green’s function was obtained using the re¬ 
cursive Green’s function algorithrrPH and the self-energies 
were calculated using the decimation methocP^. The real 
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FIG. 4. Modeling the transition from ballistic to diffusive 
transport in a micron long graphene, (a-b) sample potential 
landscape for various impurity concentrations, (c) total con¬ 
ductivity as a function of channel carrier density from ballistic 
(cr oc n) to diffusive (cr oc ^/n). 


space Hamiltonian matrix for the channel and the con¬ 
tacts were obtained using Eq. (§ and the effect of mag- 
netic field was included using the Peierls substitutiorP^. 
To calculate the voltage at contact c, we set la = —h = 

Ic = 0 and /jLa = Ey then solve Eq. 0 for lie and fiij 
where Ey is the Eermi level obtained from the electron 
density n. Then, V = l^c — and resistance R = V/I. 
The whole procedure was repeated for each n and B. 

Eig. [^c) shows the Resistance R as a function of the 
electron density n and magnetic field B. The color map 
was generated using the quantum mechanical (NEGE) 
approach and the the dashed lines were computed using 
the semi-classical formula Eq. The dashed lines and 
the bright bands represent the focusing of electrons to 
contact c. These results agree well with the experiment^. 

As a final example, we show how our model can in¬ 
terpolate between the ballistic and diffusive limits. We 
model the transport with impurity scattering by using 
a sequence of Gaussian potential profiles for the scat¬ 


tering centerP^, U{r) = Yln=i f/„exp(-|r -r„|^/2C^) 
that specifies the strength of the impurity potential at 
atomic site r. are the positions of the impurity atoms 
and ( is the screening length 3nm for long range scat- 
terers). The amplitudes Un are random numbers follow¬ 
ing a Gaussian profilP^, Wmp is the impurity concentra¬ 
tion. With U added to E[ (potential landscape shown in 
Eig. l^-b), we study the evolution of electron transport 
in graphene from ballistic to diffusive (for varying A^imp- 
Fig.|§ shows the graphene conductivity at various chan¬ 
nel carrier densities for several impurity concentrations. 
This time we keep the two ends of the graphene device at 
constant doping to capture the contact induced doping 
in graphene. This produces electron-hole asymmetry in 
the ballistic limit due to formation of pn junction near 
the contact. At high impurity concentration, the contact 
resistance becomes less dominant compared to the device 
resistance and the electron-hole asymmetry washes out, 
similar to what is seen in experimentP^. Eurthermore 
at the diffusive limit, a becomes proportional to n for 
a sample dominated by long range scatterers and can be 
fitted with a carrier density (n) independent mobility. As 
we lower the impurity concentration towards the ballis¬ 
tic limit, the graphene conductance (G) becomes propor¬ 
tional to the effective doping Ep, G = Aq^WEyT / (nfiVY) 
leading to a sub-linear a oc \/n), where the transmission 
T is determined by the metal-graphene contact. Such 
an evolution of electron transport in graphene has been 
verified in experiments^. 

In conclusion, we have shown that for massless Dirac 
fermions, an additional quadratic term in the Dirac 
Hamiltonian not only circumvents the Eermion doubling 
problem in a spatial lattice but also has a huge computa¬ 
tional advantage over the atomistic tight binding model. 
In particular, we have shown that the modified Hamilto¬ 
nian results in an extremely small matrix on a real space 
square lattice. As a result, the Hamiltonian is orders 
of magnitude faster than the tight binding Hamiltonian 
when used in band structure and quantum transport sim¬ 
ulations. We applied this Hamiltonian for micron scaled 
graphene devices to study magneto-transport and elec¬ 
tron transport in ballistic and in diffusive limit. Although 
only graphene is considered here, it is applicable to any 
other Dirac materials like topological insulators and can 
be used to calculate the spin currenfP^ as well. 
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